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Colloidal gels are out-of-equilibrium structures, made up of a rarefied network of colloidal particles. 
Comparing experiments to numerical simulations, with hydrodynamic interactions switched off, we 
demonstrate the crucial role of the solvent for gelation. Hydrodynamic interactions suppress the 
formation of larger local equilibrium structures of closed geometry, and instead lead to the formation 
of highly anisotropic threads, which promote an open gel network. We confirm these results with 
simulations which include hydrodynamics. Based on three-point correlations, we propose a scale- 
resolved quantitative measure for the anisotropy of the gel structure. We find a strong discrepancy 
for interparticle distances just under twice the particle diameter between systems with and without 
hydrodynamics, quantifying the role of hydrodynamics from a structural point of view. 

PACS numbers: 82.70.Dd; 64.70.Pf; 47.55.-t; 82.70.Gg 


The formation of a network of arrested material with 
finite zero-shear viscosity upon slight quenching is among 
the most striking features of condensed matter HHa. 
Gels can be soft or biological materials such as proteins 
[5l[6], clays [7], foods [8] , hydrogels [9] and tissues mm- 
However a more diverse range of systems including gran¬ 
ular matter m, phase-demixing oxides m and metallic 
glassformers El also exhibit gelation. The mechanical 
properties of gels are influenced by their structure both 
locally mm and at a global level through percolation 
of particles [19] and clusters [20], network topology [21] 
and confinement [22]. 

Despite its widespread occurrence, a deep understand¬ 
ing of gelation remains a challenge. In particular, numer¬ 
ical simulations with numbers of particles large enough to 
reveal the percolating network structure are limited to an 
approximation in which the effect of the fluid containing 
the particles is disregarded [23]. In the colloidal gels we 
will consider, demixing of the particles into a (colloidal) 
“gas” and “liquid” occurs. Spinodal demixing leads to a 
network of particles [24]-[28] which undergoes dynamical 
arrest [29] . The final structure can persist for years m, 
if the self-generated or gravitational stress is weaker than 
the yield stress [8] . 

Demixing is driven by effective attractions between 
the colloidal particles induced by the addition of non¬ 
absorbing polymer. Thus although the system is a mix¬ 
ture of three important components — colloids, polymers 
and solvent — in equilibrium we can recast it as an ef¬ 
fective one-component system of colloids which experi¬ 
ence an attractive interaction whose strength is deter¬ 
mined by the polymer concentration [SUEg. However 
out of equilibrium, hydrodynamic interactions (HI), i.e. 
interactions between particles mediated by the continu¬ 
ous fluid phase, come into play. These are hard to treat 


as they are both long-ranged and not limited to pairwise 
interactions. 

In this Letter, we show that the colloid volume frac¬ 
tion required for gelation is reduced by a factor of 1.75 
including HI, see Fig. [^ This occurs although our Brow¬ 
nian dynamics simulations (without HI) are parameter¬ 
ized closely to our experiments [33l [34] . We back up these 
findings with simulations including HI using fluid particle 
dynamics [35l [36] for accessible system sizes. Our results 
call into question the standard theoretical approach used 
to describe colloidal systems quantitatively. 

We show that HI strongly favor gelation since they 
inhibit the direct formation of energetically locally fa¬ 
vored compact structures, which minimize the potential 
energy between small groups of m colloids [371138]. On 
one hand, the structural mechanism by which the net¬ 
work of colloids forms and becomes rigid has been shown 
to be condensation into locally favored structures whose 
geometry underlies the local order of the percolating net¬ 
work [34]. On the other hand, HI lead to solvent back 
flow, with the effect of slowing down the condensation 
of colloids, favoring more open elongated local structures 
than would be the case without HI (see Fig. inset) 
[351 IMi I 39 H 4 T] . This is a consequence of the incompress¬ 
ibility of the liquid solvent, which allows only transverse 
(rotational) motion of the fluid. Without HI, large closed 
structures form, favored by minimizing the energy alone. 

We used polymethylmethacrylate colloids, sterically 
stabilized with poly hydroxyl steric acid. The colloids 
were labeled with 4-chloro-7-nitrobenzo-2-oxa-l,3-diazol 
(NBD) and had a diameter of a = 2.40 jam and poly- 
dispersity of 4% as determined by static light scattering 
[42] , The polymer was polystyrene, with a molecular 
weight of 3.1 X 10^ with M^/Mn = 1.3 and solvent was a 
nearly density- and refractive index-matching mixture of 
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FIG. 1: (color online) State diagram showing gelation in our 
experimental and simulated system in the f3e (potential well 
depth)-0 (colloid volume fraction) plane. Phase separation 
is observed for /3e > 3. Filled symbols are experimental 
data, Brownian dynamics simulation data are unfilled sym¬ 
bols. Red symbols are one-phase fluids, circles lie below the 
percolation threshold, triangles above. Black circles are (non¬ 
percolating) cluster fluids. Blue squares are gels. Inset shows 
effects of hydrodynamics on colloidal aggregation [39]. Top 
line: hydrodynamic interactions lead to solvent flow (cyan 
arrows) which influences the motion of the aggregating col¬ 
loids. The incompressibility of the fluid allows only transverse 
(rotational) flow, resulting in the formation of an elongated 
structure rather than a closed one. The degree of elongation 
is characterized by bond angle 6. Bottom line: no hydro- 
dynamic interactions. Their attractive interactions lead the 
colloids to aggregate (pink arrows) to form a compact struc¬ 
ture. 


cis-decalin and cyclohexyl bromide. To screen any (weak) 
electrostatic interactions, we dissolved tetra-butyl ammo¬ 
nium bromide (TBAB) salt to a concentration of 4 mM 
[33l [34] . We determine the colloid volume fraction 0 and 
polymer number density pp by weighing along with par¬ 
ticle tracking in the case of the former. 

We take the polymer volume fraction as <pp = 
A'KR^gppl3. In these good solvent conditions Rg ^ 160 
nm so the polymer-colloid size ratio is g ^ 0.18 [33] . This 
is much larger than the range of the steric stabilisation 
layer (^ 10 nm) previously identified as a criterion for 
spherically symmetric effective interactions between col¬ 
loids, which we assume here, neglecting for example any 
hindrance of rotational diffusion [43] which might other¬ 
wise impact significantly on relaxation dynamics [44] . 

We prepare the samples in glass capillaries of inter¬ 
nal dimensions 100 pm x 1 mm x 1 cm and image the 
colloidal particles with a Leica SP5 confocal microscope 
with a resonant scanner and took data at least 20 pm 
from the wall of the sample cell. The particle coordinates 
are tracked in three dimensions (3D) with an accuracy of 


around 0.03cr. Samples were prepared 10 minutes prior 
to imaging and imaged for around 20 minutes. In units 
of the Brownian time, tb = {al2)‘^ jdD = 3.1 s, where 
D is the diffusion constant. Data was therefore taken 
between 192 and blSrB rafter preparation. In our simu¬ 
lations we model the colloid-polymer mixture, using the 
Morse potential 

0u{r) = ^eexp [pQ{a - r)] (exp [po{a - r)] - 2) (1) 

which we have found to provide an accurate description 
of colloid-polymer mixtures for size ratios such as we con¬ 
sider here [34l [45] . 

The interaction is truncated and shifted at r = 1.4cr 
and the range parameter is set to po = 33.0 [34]. The 
effective temperature is controlled by varying the well 
depth of the potential scaled by the thermal energy 
k^T = 1//3, Pe, which corresponds to the polymer chem¬ 
ical potential in the experiments. This can be estimated 
from the polymer concentration using Widom particle 
insertion imiiT]. Short-ranged attractive systems can 
be mapped onto one another [48l |49] via the extended 
law of corresponding states, which equates the reduced 
second virial coefficient = 2B2/{37ra^f^) between two 
potentials where aeff is the effective hard sphere diame¬ 
ter EH Ellis]- Here B 2 is the second virial coefficient. 
Knowing the polymer volume fraction in the experiment, 
we are able to map the effective colloid-colloid interaction 
to a Morse potential of the same second virial coefficient. 
We assume that the colloid-colloid interaction has the 
so-called Asakura-Oosawa form [33]. Further details are 
provided in Supplementary Material (SM). 

Our colloidal particles used in the experiments are 
slightly polydisperse, which is treated in our Brownian 
Dynamics (BD) simulations by scaling r in Eq. 0 by 
a Gaussian distribution in a with 4% standard deviation 
(the same value as the experimental system). For details 
of our BD simulation methods, see [50]. Unless otherwise 
indicated, runs are equilibrated for WStb (close to the 
value in the experiments). Data are then sampled for 
a further WStb- The system size was fixed by the box 
length of 47.13 a up to = 40000 particles (correspond¬ 
ing to (j) = 0.2). For higher volume fractions, the box 
was reduced and the system size fixed at A" = 40000. In 
the simulations we define (j) = Acr^7r/(67r). To analyze 
the local structure, we identify the bond network using 
the Voronoi construction with a maximum bond length 
of 1.4(7. Having identified the bond network, we use the 
Topological Cluster Classification (TCC) [38] to decom¬ 
pose the system into a set of locally favored structures 
comprised of m particles, which are the minimum energy 
clusters for the Morse potential and illustrated in 
Fig. [^a) and (b) [38] , 

The most striking observation in the state diagram of 
Fig. 0 is that experiments form gels (in the sense of 
a percolating network) even at 0 = 0.04 ± 0.005 [51] . 
By contrast BD simulations form gels only above (j) = 
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FIG. 2: (color online) The effect of hydrodynamics on the local structure, (a-b) Snapshots illustrating the effect of hydrody¬ 
namics on the local structure of non-equilibrium sticky spheres. Two nearly identical state points are shown, with particles 
colored according to their local environment as identified with the TCC for 0 = 0.05, (a) Experiment, fie — 8.3. (b) BD 
simulation, [Se = 8.0. Cluster topologies are indicated at the right of (b). (c-d) The effect of hydrodynamics on the local 
structure for 0 = 0.05. At low interaction strengths there is little clustering, (c) Experimental data are dominated by m = 5 
triangular bipyramids with small quantities of higher-order clusters, (d) BD simulation data show the assembly of many more 
clusters of more complex geometries. Vertical dashed lines in (c,d) correspond to the (critical) attraction strength required for 
demixing and line labels refer to the number of particles in the cluster. 
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EIG. 3: (color online) Population of clusters identified with 
the topological cluster classification. Experimental data cor¬ 
responds to 13s = 8.3, for simulations jSs = 8.0. Simulation 
data is shown for fluid particle dynamics (HI) for N = 6112, 
Brownian dynamics for a larger polydisperse system (BD) 
and a smaller monodisperse BD system with N — 6112 (BD 
mono). Highlighted are m = 4 tetrahedra and m = 5 triangu¬ 
lar bipyramids favoured in the case of HI and m = 10 clusters 
favored by BD simulations. 


0.07 ± 0.005, requiring a much higher volume fraction 
than in the experiments. We consider different simulation 
run times and system sizes, neither of which have any 
significant effect on our results (see SM). 

We use the TCC to identify different locally favored 
structures and consider a state point which typifies the 
effects of HI on local structure in Fig. At ~ 8, 
the experiment exhibits a percolating gel, while BD sim¬ 
ulations form isolated (non-percolating) clusters. Un¬ 
der these conditions, the difference in local structure 
is clear: the experiment (a) is dominated by smaller 
clusters, mainly tetrahedra and five-membered triangu¬ 
lar bipyramids. By contrast the simulation (b) exhibits 
larger clusters. In Figs. |^c) (experiment) and (d) (sim¬ 
ulation) we consider a range of well depths. We see that 
over the range relevant to gelation, around 35% of the 
system is in tetrahedra in the experiment, with a fur¬ 
ther 35% is in m = 5 triangular bipyramids, which are 


formed by the addition of a particle to the center of a face 
of a tetrahedron; these two simple clusters dominate the 
experimental system. This contrasts strongly with the 
BD simulation data, where much larger, compact clus¬ 
ters based on five-membered rings (m = 8, 9 and 10) 
account for over 60% of the particles in the system. Note 
that tetrahedra are the smallest rigid unit for spheres in 
3D, thus we expect that these represent a limiting case 
for a dynamically arrested gel structure. Our findings for 
the local structure are thus consistent with the schematic 
picture outlined above and also with previous simulation 
work which showed that HI can suppress the formation 
of compact clusters [35l [36l I39H4T] . 

To provide further evidence in support of our hypoth¬ 
esis we carry out fluid particle dynamics simulations 
(FPD). This method includes HI by treating both col¬ 
loids and solvent as a fluid continuum with a sufficient 
viscosity contrast [35j |36], to realize a situation where 
viscous dissipation occurs predominantly in the solvent 
as appropriate to a suspension of solid colloidal parti¬ 
cles. Such simulations necessarily describe smaller sys¬ 
tems due to the high computational cost. Here we con¬ 
sider a monodisperse system of N = 6112 particles at 
a volume fraction of (j) = 0.05 and attractive well depth 
/3s = 8. The results are shown in Fig. and Fig. 
Hd). We also carried out a monodisperse BD simulation 
with N = 6112, i.e. matched to the FPD parameters. 
We see that FPD simulations with HI are much closer to 
the experimental results than either BD result, indicating 
that the main difference between the BD simulations and 
experiments is HI, rather than polydispersity or system 
size. 

Thus we argue that the origin of the percolation at low 
colloid density in the experiment is the lack of compact 
clusters which lead to a more open structure. These more 
open clusters arise from hydrodynamic back flow during 



































4 


condensation. While the FPD simulations are very much 
closer to the experiments than the BD simulations, the 
agreement is not perfect. This might originate from dif¬ 
ferences in the system size, as the experimental system is 
very much larger. Furthermore the boundary conditions 
imposed by the capillaries in the experiment are rather 
different to the case of of the periodic simulation box. 



FIG. 4: (color online) Three-body correlation functions to 
probe hydrodynamic effects gsk as defined in the text; 0 is 
the bond angle, (a) Experiment for /3e = 8.3. Inset: FPD 
simulation for /3e = 8.0. (b) BD simulation for f3e = 8.0. 

Inset shows the same data expanded to show the full extent 
of the peak around \/2cr. Simulation data has Gaussian noise 
(width O.lcr) added to the coordinates to mimic experimental 
errors m- 

Having identified the fundamental local structures, we 
now investigate the lengthscales on which HI lead to 
thread-like structures, which make up the gel network. 
Inspired by similar efforts to detect elongated structures 
in cosmology [52] , we consider three-body correlations as 
function of the bond angle. Figure[^ (inset) indicates that 
the bond angle distribution is likely to be biased toward 
bond angles 0 close to tt compared to the case without 
HI |36l [41] . A general three-body correlation function for 
a triad r, r', r — r' may be written as [53] 

fl 3 (r,r') = 6{rk - Vi + r)6{rk - rj + r')) , 

P i^j^k 

where p is the number density, and ijk are particle in¬ 
dices. By considering g^{r^0)^ where |r| = |r'| = r and 
0 is the angle between |r| and |r'|, we probe the angle 
at different length scales r. For simplicity, we compare 
0 = IT and 7r/2, a strong signal in the former demonstrat¬ 
ing elongated structures typical of a rarefied gel. 

We remove the contribution of pair correlations, which 
can mask the differences in bond angles we seek, using 
Kirkwood’s approximation ^ 3 (r,r') ^ ^(|r|)^(|r'|)^(|r — 
r'l). We thus plot gsKir,!^) = 9s{r,T^)/[g{r)g{r)g{2r)] 
and 5i3K(r,7r/2) = g:i{r,'Kl2)/[{g{r)g{r)g{^r)\ at (5e « 
8 for both experiments and simulations in Figs. (a) 
and (h) respectively. We see in the simulation data in 
Fig. mb) that there is a very strong peak at r ~ \/2cr for 
0 = 7r/2, corresponding to the compact structure shown 


in the inset. As expected, no corresponding signal is seen 
at 6> = TT, which is sensitive to elongated structures. 

By contrast, there is a complete absence of any peak in 
the 0 = 7t/2 data of the experiment, as seen in Fig.j^a). 
Such a reduction in compact configurations indicates that 
the presence of HI has a dramatic effect on the geometry 
of configurations at the three-body level. On the other 
hand, experiments show a significant peak at around r ^ 
1.7a < 2(7 for 0 = it, [cf. Fig. |^a)]. This is indicative 
of chainlike structures with a length scale ~ 3 — 4(7, as 
expected in the case of HI, which tend to align colloids 
along the flow, as shown schematically in Fig. [^ There 
is no detectable signal from clusters such as defective 
icosahedra, whose length scale is ^ 3 — 4(7 as well, as 
they do not produce chains sufficiently straight to lead 
to a peak for 0 = tt. As shown in Fig. I^a) inset, FPD 
simulations closely follow experiments in their three-body 
correlations. 

In conclusion, we have shown that hydrodynamic in¬ 
teractions have a profound effect on gelation in colloids: 
the volume fraction above which gelation is found is 
(j) ^ 0.07 for simulations without hydrodynamic inter¬ 
actions whereas cj) ~ 0.04 for experiments with HI. Fur¬ 
thermore, we have found that the effects of HI can also be 
seen from the distinct difference in the particle-level local 
structures formed upon phase demixing between experi¬ 
ments and simulations without HI, which we have backed 
up with simulations which include HI. In the absence of 
HI, compact clusters comprising around ten particles are 
formed. However those formed in the case of HI at the 
same volume fraction and attractive strength are pre¬ 
dominantly four-membered tetrahedra or five-membered 
triangular bipyramids. Tetrahedra are the minimum me¬ 
chanically stable structure. Thus HI promote the for¬ 
mation of tenuous low-density networks. The range over 
which HI impact local structure is determined with a 
three-body correlation function. This is sensitive to the 
spatial extent of anisotropic bond angle corrections. We 
find that the effects of HI can be observed over distances 
of 4(7 or around 10 pm at 0 = 0.05. We show that, when 
assessing the low-density limits of gelation, it is crucial 
to consider hydrodynamics. Thus in designing products 
based on gels where the concentration of gelling agent is 
important, one should carry out experiments or simula¬ 
tions which include HI, rather than to rely on Brownian 
dynamics simulations. 
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Supplementary Material 

Mapping interactions between experiment and simula¬ 
tion — . We have previously [33] demonstrated that this 
system is accurately described by the Asakura-Oosawa 
interaction between the colloids, which reads 


pUAo{r) = < 




1 - 


3r 


2(l+<3f)cr ' 2(l+g)3cr3 


0 


for r < (T 

for (7 < r < (T ^ (Tp 

for cr -h cTp < r 


( 2 ) 



FIG. S5: (color online) Reproducing experiments with simu¬ 
lation. The pair correlation function in experiment (circles) 
is reproduced by Monte Carlo simulations with the Morse po¬ 
tential. These have a polydisperse colloid size distribution of 
the same form as the Brownian Dynamics simulations. In the 
simulations the well depth of the Morse potential /3e = 1.8. 

Here (jy^ is the polymer volume fraction in a reservoir in 
chemical equilibrium with the experiment. Widom par¬ 
ticle insertion @ 61117 ] is used to map between reservoir 
and the experimental system, cjp = 2Rg is the poly¬ 
mer diameter where Rg is the radius of gyration. The 
Asakura-Oosawa interaction is mapped to the Morse po¬ 
tential as described in the main text. 

To investigate the accuracy of our technique in Fig. 
we show the pair correlation function g{r) of a stable fluid 
reproduced with simulation. A Gaussian distribution of 
standard deviation of 0.03cr is added to the coordinates to 
mimic the particle tracking errors in the experiment [33] . 
The agreement in the g{r) indicates that the pairwise 
interactions in the experiment and simulation are very 
similar. 
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other symbols are the same as Fig. (1) in the main text. 
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